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ABSTRACT 


A UNIVAC 1108 Computer has been used to obtain exact values for _ 


{ip 1, 2, ..., N), the non-zero characteristic values of Zimm's matrix HA, 


for values of N (the number of Gaussian springs per molecule) up to 300; 


valies previously published did not go above N = 15 (Thurston and Morrison), 


fo future comparison with elastic and viscous properties of very dilute polymer 


solutions measured in steady and in ¢cscillatory shear flow, several useful 


functions of the \_ are computed and tabulated for parameter values in the 


1/2 


ranges 2<N< 300, 0.05<h <0.3, where h (=hN~’“%) is a dimensionless 


"friction constant"; h denotes the "hydrodynamic interaction parameter” of 


Zimm, For selected parameter values, certain of the functions are compared 


graphically with the corresponding functions given by Tschoegl on the basis of 


approximate calculations using an indefinitely large value for N. The present 


calculations differ from those of Zimm and Tschoegl in only two ways: finite 


values of N are used (instead of N = ©), and exact values of Me are 


computed irom an N XN _ matrix equation (instead of approximate values from 


an approximate integro-differential equation). An important new result of the 


present computation is that, for the particular value h - 0.26 2, Plorv's 


intrinsic viscosity function ® is substantially independent of N in the range 


3 
50 < N < 300 and there has the value 2.85 X 107 ml/g; this equals the value 


given by Zimm's original approximate calculation forthe case N “™, h>>l. : 


This suggests that the interpretation of intrinsic viscosity data in terms of the 


well-known ‘non-tree-draining' concept should be reconsidered, 
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EXACT RELAXATION TIMES AND DYNAMIC FUNCTIONS FOR DILUTE POLYMER 
SOLUTIONS FROM THE BEAD/SPRING MODEL OF ROUSE AND ZIMM 


i A. S. Lodge and Yeen-Jing Wu 


4 1. Introduction 
In the well-known theories of Rouse (1) and Zimm (2) for very dilute 
solutions of long-chain (‘linear polymer’) molecules, each polymer molecule 
is represented by N equal ‘therma! motion’ ‘Gaussian’ springs connected 
end-to-end; the solvent is treated as a constant-temperature environment for 
the polymer molecules and also as a Newtonian incompressible liquid, of 
viscosity Ns which interacts with each polymer molecule only at N +1 point 
centers ('beads') situated at the ends of the N springs. Lach bead has the 
, same ‘friction constant’ fy (which would be equel to 6ran. if the bead were 
a rigid sphere of radius a treated according tu Stokes‘ law). Each spring 
has the same root-mean-square end-to-end distance hb when the solution is 
not flowing. 
Compared with earlier theories (3) 1n which equal friction constants were 
‘ assigned to each atom on the ‘backbone’ of a polymer molecule, the bead/spring 
theory of Rouse and Zimm involves an extra artificial feature, namely, the 
concentration of polymer-solvent hydrodynamic interaction at the N+ ] 
fictitious beads. The advantages arising from the simpler mathematics involved, 


however, are considerable: the dynamic viscosity n' and dynamic rigidity G' 
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* Presented at the Society of Rheology's Winter Meeting, Salt Lake City, 
February, 1971. 
Sponsored by the United States Army under Contract No. DA-31-124-ARO-D-462., 
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as functions of angular frequency ~, for oscillatory small-amplitude shear, 
have been calculated (1, 2, 4) and compared rather extensively with experi- 
mental data (5); more recently, the (macroscopic) constitutive equations, from 
which all rheological properties in isothermal homogeneous flows can be 
calculated, have been derived by Lodge and Wu (6), who also list all the 
assumptions involved in the theory. Precise calculation of the dynamic functions 
n‘(w), G'(w! has recently assumed added importance because improvements in 
experimental techniques have made possible the extrapolation to zero polymer 
concentration of the data obtained (7, 8); such ‘intrinsic’ data can thus properly 
be compared with the predictions of the bead/spring theory which is based on 
neglect of interactions between different polymer molecules. 


In order to compute the functions n'(w), G'(w), it is necessary to compute 


the characteristic values - (p= 1, 2, ..., N) of the symmetric matrix N XN ; 
B { Bod! whose elements Bog are given by the equations 
7H +H - H -H S41 By ead ol 
Boa pq p-l,q-l p-l,q p.q-l (P, » Gy , N), {1 J 
where 
1 (i =j) ; 
Hay = 4 oa i (ij $30, 2s sce ND, [1.2] ; 
an lies |? tig) : 
; 
Ba 3 1/2 ' j 
ho = fp/[U2m yp" nob] . [1.3] 
i 
j 
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The matrix B has been introduced by Iodge and Wu (6) as a more convenient 
alternative to Zimm's singular, non-symmetric, matrix HA which has the 


characteristic values 0, \ Ny acer 


N (6); H “LHI, and A is defined 


in equation (ll) of Zimm (2). 
The object of the present paper is to present the results of computing 


‘exact’ characteristic values Aya Non vey ® of the matrix B, tcgether with 


N 
various functions of the - of interest for comparison with experimental data 
tained in steady- and oscillatory-shear flow. Available computer storage 


facilities have enabled us to use values of N_ up to 300, which may well prove 


to be sufficient for the highest-molecuiar-weight polymers available. The only 


other published exact computations (for non-zero h ) are those of Thurston and 


Morrison (9), whose restriction N < 15 is rather severe. 
We call our values for Me ‘exact' to distinguish between them and the 


approximate values obtained by Zimm, Roe, and Epstein (10), Hearst (11), and 


Tschoegl (4) (Table 1), who replaced the discrete eigenvalue preblem (associated 


with matrix equations involving B or HA) by a discrete eigenvalue problem 
associated with an integro-ditferential equation (equation (71) of Zimm (2)) 
regarded as an approximation to the actual matrix equation given by the Zimm 


model, The approximation involved here has becn thought to be reliable when 


N is large; recently, however, partly as a result of comparing our exact values 


for aS with those obtained from the Zimm integro-differential equation, 


Osaki (12) has noted that Zimm's equation is incorrectly derived from the matrix 


equation, and that the corrected equation is the following: 
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ie 1 
a(r(l-4h th f a(s)|r- s| > -Nda(r)/4 . [1.4] 
zi 


* 
Zimm's equation omits the term -4h a(r). Zimm's parameter h is given by 


the equation 


neh’n’? , [1.5] 


In Table 1, the notation h >>1 (used by Zimm) mears that, on the left-hand 


: side of [1.4], the integral term alone is retained; as Osaki (12) has observed, 


this case really corresponds to the value h- - 0.25 (and N large) when the 
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correct approximate equation [1.4] is used. 
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q TABLE |: Calculations of characteristic values _ for the 
7 bead/spring theory of Rouse and Zimm 5 
q N h h” ; 
a 1, Rouse (1) 1 to © 0 0 exact : 
4 : 
E: 2. Zimm, Roe, and Epstein (10) large >>1 approximate i 
q 3. Hearst (11) large >0 approximate : 
1 
z. 4. Tschoeg! (4) large >0 approximate é 
a 5. Thurston and Morrison (9) 1to15 0 to 0.4 exact : 
i Ff 
E 6. Wu (13) lto 222 0.05to0.3 exact 5 
7. Lodge and Wu (this paper) 2 to 300 0.05 to 0.3 exact : 
3 ; 
i In view of the idealization inherent in the bead/spring mode] and 3 
; various approximations used in Zimm's analysis of this model (on which our 
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present calculations and all the calculations of Table | are based), it might 
be argued that the present computations amount to hair splitting; we recognize 
the weaknesses that remain, but consider that it is worth while making the present 
computations if only to eliminate any easily avoidable approximations in the 
analysis of the model and so render it more open to critical examination. Pyun and 
fixman (14) have made certain calculations based on the bead/spring model 
without using one of Zimm's assumptions (namely, that the Oseen hydrodynamic 
interaction tensor can be replaced by its value averaged in the solution a* rest); 
Wu (13) has compared results of Pyun and Fixman's calculations with the results 
of the present paper. 

The characteristic values te are dimensionless functions of the 
parameters N, h. Certain constants To called 'relaxation times', defined 


; by the equations 


= (m 73/2 b*h’/(kT) (p =1, 2, «+e, N), [1.6] 


where k and T denote Boltzmann's constant and absolute temperature, play 
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a fundamental role in the theory. The present report thus computes exact i 


amet 
tay 


relaxation times in terms of three parameters, whose values are unknown 
a priori, The selection of these three parameters may be made in various ways, ; 


some of which are listed in Table 2. 
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TABLE 2. Selections of three basic parameters 


ly. N,f., ) N_ springs per molecule; fy = bead friction constant; 


b =rem.s. spring length at rest; 


x * 1 
Ze. (Nhe gb he fy/| (an?) 2, bhi 
3, Njh,7, 7, 2 fb°/L6 TAN, 1; 
* 
4, N,h,7, nen'nY 4, 
& N 
5. N,h,7 T, = Y Tes 
’ 7 "Q 0 = p 
‘ = 
| 6. m,h,b m= M/N (M = polymer molecular weight). 


We follow the usual notation according to which \ i* 5 Sele < Ay ° 


\, is thus the smallest characteristic value, and hence 7, is the greatest 


1 


relaxation time. 


2. Computation of Np ros eoey An ° 


i Ay Aas eees My are the roots in X of the equation 
i 
oy det (B._ - 46) =0. [3.1] : 
: pq pg 4 
4 When h- = 0, the roots are given by the equations ‘ 
oad sin@| pr/2(N + 1)] (p=1, 2, «+, N) [3.2] 


*% 
given by Rouse (1). When h #0, the roots for the cases N =1, 2 areas 


a follows: 
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dped-b, 42340 2/24yn (N= 2). 


For h #0 and N> 3, we know of no method of obtaining a solution 
for the \_ in closed form. To obtain solutions by numerical methods, a Fortran 
program was written in order to compute roots of [ 3.1] by means of a combination 
of Householder's method and the 'Q-R method'; the University of Wisconsin- 
Madison Computing Center's UNIVAC 1108 and ‘BUMP 2°‘ Library Subroutine were 
used; in this way, values of Ap ho» exes aN were obtained for various values 
of i and for values of N up to 222 (13); computer storage limitations made 
it impossible to handle higher values of N using this program. Subsequently, 
use was made of the fact that B is symmetric, and a modified Miami University 
Program (H. J. Wertz version, dated 4/5/68) for handling packed matrices 
enabled the computations to be extended up to values of N = 300. To compute 


a complete set of characteristic values for a given value of N anda given 


value of h , the following times were required: 


N: 25 50 100 « O 300 


Time in seconds: 0.3 1.6 10 70 235 . 


* 

These times do not vary much with h . 
Table 3 contains the results of computing va~ious functions of the x ‘ 
The values of A, and dh. are tabulated. Values Of day day very Nyy have 


not been tabulated here, but are available on request to one of the authors 


(A. S. L.); some have been tabulated by Wu (13). 
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The following checks on the computation confirm the reliabilit. of the 
results and indicate that the computational error in Table 3 entries should 
nowhere exceed unity in the fourth significant figure: 


1, Values obtained for NS are consistent with the equation 


= Trace B; L357} 


Timez 
” 


l p 

2. For N =1, 2, ..., 15, and various values o1 h, the values 
obtained for ‘ agree with those computed by Thurston and Morrisun (9); 

3. For N=l1, 2, values of Ms are consistent with [ 3.3]; 

4. For N = 50, 100, 200 and h . 0.1, 0.2, 0.3, most of the 
tabulated functions were computed using both progranis referred to above; 
comparison of entries showed agreement to three significant figures and, in 
most cases examined, to four significant figures; a few entries showed 
differences of unity in the fourth significant figure. 

A minor advantage in using h’, rather than h, as a parameter in 
the tabulation is that h- has an upper limit, while h does not (assuming 
that N could be indefinitely great). We have chosen values of h” in the 
interval 0.05 < h* < 0.3. Thurston and Morrison (9) quote the value 0. 471 
for the upper bound of h’, but do not give details of their derivation; in an 
attempt to follow their argument, we obtain the value 0.65. Alternatively, 
using the very crude argument that fy = 6man . for a rigid spherical bead of 


radius a, the inequality 2a <b yields the result h <0.49. Because of 
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the artificiality inherent in the model, we do not piace much weight on the 

: actual value of the upper bound on bh” given by such arguments; it seems 
reasonable, however, to suppose that there is an upper bound. Independent 
arguments based on a compariscn of the results of the present calculations 
with measured data on intrinsic viscosity of polymer solutions suggest that 
the useful range for hg is 0< h < 0.26. This range is fairly well covered 

by the entries in Table 3. Thurstor and Morrison (9) have noted that the 


condition h” < 0.65 is sufficient to ensure that a shall be positive for 


N-1, 2 (cf. [3.3]). 


3. Steady shear flow functions 


In the usual notation, we express the intrinsic viscosity [yn] in the form 


n-7n. : 3 set 
8 @ cpty3/2y ] (aii ** 


[yn] := Jim 
c-0 gs 


where c is the polymer concentration in g/ml, n is the solution viscosity, 
and <r°> is the mean square end-to-end distance of a polymer molecule in 
the soiution at rest. This equation may be regarued as a definition of the 
quantity ®, ‘The bead/spring theory gives the following equation for 
® = O(N,h): 
as (3/52 nw 3/2 9 s [3.2] 
p=l 


where N > 6.0228 X 1073 (Avogadro's number). 


ate ate 
m8 5,8 


The notation A :=B rreans "A is defined by the equation A = B", 
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In Table 3, values of ® are tabulated for various values of N and h. 
for comparison with the approximate values calculated by Tschoegl for N =, 
the values of © have been plotted as functions of logh in Figure l. It 
is seen that the experimentally - inte: 2sting range of values for ®, namely 


2.5 %10 °° 


@ < 2.9, are obtained with values of N inthe range 25<N < 300 
and values of h near 6.25, when the exact characteristic values are used, 
and with large values of h when the approximate characteristic values (with 


N =) are used. In particular, for the special value h- = 0,262, Table 3 


gives the following values for ® (in g/ml): 


N: 8 25 50 100 200 300 


[3.3] 
10°°7@: 2.954 2.864 2.853 2.852 2.855 2.856 (h =0.262) 


For 25<N < 300, it is seenthat ®(N, h) (with h’ : 0.262) is substantially 
independent of N and that its value is close to the value 2.84 x 1073 obtained 
by Zimm (2) from the approximate calculation based on the values N = ©, 
h>>1l, 

The case h >>1 is described by Zimm (2) as the "non-free-draining 
case", presumably because a large value of h is associated with a predominating 
bead-bead hydrodynamic interaction term in the equations of motion for a bead/ 
spring unit and this in turn is believed to imply that, within a surface enveloping 
the polymer molecule, there is a substantial region of solvent which is 


‘immobilized’ by the polymer molecule in the sense that the solvent velocity 
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equals some average polymer molecule velocity. Whether or not this concept 
could be substantiated tor the case h >>} by calculation of selvent and polymer 
velocities using the equations of the Zimm theory, our present results suggest 
that the concept of the non-free-draining limit may not in fact be appropriate 
for dilute solutions of linear polymers for which ©@ has a value near 
2.85 X 1073 g/ml: using our exact characteristic values, such a value is 
obtained when h” - 0.262, and there seems to be no a priori reason for thinking 
that this value of the dimensionless friction constant is so iarge that partial or 
complete ‘solvent immobilization’ should occur. 

In Figure 2, the data are replotted to show the variation of log ® with 
h’ at constant N; the fact that the curves for N > 25 pass approximately 
through a common point (whose abscissa is 0. 262) reflects the results [ 3.3]. 

For comparison of the theory with the experimentally determined molecular 
weight dependence of the intrinsic viscosity, Thurston und Morrison (9) make 
the reasonable assumption that N is proportional to M. For limited ranges of 
M values, both [yn] and <r’> (determined from light-scattering data) vary 
epproximately as powers of M. In this connection, it is of interest to note that, 
according to our calculations, ® can be represented approximately by an equation 
of the form 


. nh) 


®(N, h) = a,(h)N : [3.4] 


for limited ranges of values o1 N; the dependence of the index ny on h 


is sho.vn in Figure 3: the fact that two curves are given for different ranges of 
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Vigure 2. Variation of log ® with h at constant N calculated 


from | 3.2] using exact characteristic values Ns . 
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N reflects the fact that | 3.4; 1s only an approximate representation of the 


dependence of ® on N. Again, the fact that n, = 0 when h - 0.262 isa 


1 


reflection of { 3.3]. 


In Figures 4, 5, the quantities 


V \ 
tft um Te = Ls d/o [3.4] 
p= p= 
N N 
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J _pel - _ pal [3.6] 
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are represented as functions of log h at constant values of N, and are compared 
with the corresponding quantities evaluated by Tschoegl (for N . %) using 


approximate characteristic values. 
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Figure 3. Dependence of index n, on h , where ®(N,h ), 


1 


calculated from {3.2} using exact characteristic 


values p? is represented approximately by [ 3.4]. 


S162 


ggg es eT ee a Ae a oe 


oo Is ; 


(N/,4 = YU) 


*(o = N) [beoyosy 


jo suoyzernoyeos ayewtxosdde yy (002 ‘00T 


‘G2 = N) SUOTJETNOTeS JORxe Jo uOSTIedWIOD ‘N 


ted 
WUEYSUOD ye Y GOT YYM TAs jo uopyeyien 96°“ SaNbTg 
N 


(|B20YDSL) 00 


een 


#1250 


} 
seg ia Costas 2 cried 
als dt. ARE We AE, G 


ALEC TORN T StS fet AD ADS) Maden iiroeere * — ~ 
A " ha alot an Th OMS «5 yer 


Fetes 


— = 6 g ‘ a 
SPIES ACRE I NEA TREY LEP EE NES IO NTT LOLA TI NE ie UNITern tiene me ss 7 


a 
Ls} 
ZO 
! oz 
. ° (o = N) ,baoyosy \ ay 
3 c 
: Jo uoyefnoyeo ayeuyxoidde ym (002 ‘0S “OT = N) we dy a 
q SUOFZETNOTeO exe jo rie egal *N jUeYSUOD . o 
E ye y 60 A (Sa ) youopeyen “SG smbTT so. ; 
| y So, wt Obs K ee _v0 
P| 
¥ ; . m4 
a = w 
E St) 
q se 
§ 7 7 vane eee oid . Te we ry isis: - Bees ~ = ~ see 


4. Small-strain oscillatory shear 


According to the bead/spring theory (2), the dynamic viscosity 1' 
and dynamic rigidity G' ina state of steady sinusoidal oscillatory shear 


of small strain amplitude satisfy the equations 


n-n & 

: S RT al 

4 lim ra, a Tmp(N, h, ow), { 4.1] 

c+0 

fi 

b G' RT “ * 

“ lim co 7M 7impN, bor), [ 4. 2] 

4 c--0 

| where R= kKN,, 

‘ | NY, 2)\-1 

: Oe ae ie) os fae 

p=l ob p 

7 Rr ut 2, [To] \b+ wT ‘ [4.4] 

2 p=l\ p p 

a and w denotes the angular frequency. The greatest relaxation time Ty 

a according to [1.6], is given by the equation 

iq 3 n 3, * 

, T~ b 

: 7) =7,(N, ae b) = )/? = oA [4.5] 

: KIN, bh) 

4 and is thus a function of the three parameters N, h, b. 

% Using the exact characteristic values, the functions TRIN ; h™ » wT ) 
4 and na(N, h, wT») have been computed for values of N, h , and log, #1 é 
i in ranges of practical interest. The results, together with values of related ; 
3 functions of interest, are given in Table 3. 
y In figures 6-9, log Rp and log TR are plotted as functions of . : 
* ‘ t 3 
4 log wr 0? where T) is defined by the equation 4 
. 4 
| 
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T =) 7 . { 4.6} 


Ty 1s a convenient dimensionless variable to choose for the abscissa when 
comparing dynamic functions TR and NR obtained using different values of 
parameters N, h- and different methods of calculation (exact and approximate): 


from [4.1], [4.3] and the fact that n =n' when w= 0, it follows that 


[n] eee [4.7] 
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Comparison of curves using the abscissa log wr, thus amounts to comparison 


0 


at constant [yn]. 

From Figure 6, it is seen that there is a small but significant difference 
between the exact and the approximate curves when h = 2; the differences at 
the largest values of log wT 1 doubtless reflect the difference in values of 
N (100, ©) used in the comparison (see also Figure 7); the difference in curves 
at lower frequencies is attributed to errors arising in the approximate calculation 


(based on the use of [1.4] without the term - 4h’ a(r)) ° 


% F 
From Figure 7, it is seen that, when h = 0.3, a change in the value of N 


BIRASY Boreseetogewenee + 


from 50 to 200 leaves the low-frequency behavior substantially unaltered but 


j 
i 


moves the point of intersection of the Np curve and the nj J-curve to higher 
frequencies. 
Figures 8 and 9 compare the Zimm curves (N = ©, h >>1) with the exact 


curves for high (h =0.3) and low (h =0.1) values of the friction constant. 
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It is intended, in a future publication, to present results of a 


comparison of the present calculations with experimental data giving 


the variation of [yn] with M andof n' and G' with wo. 


5. Discussion 

The present calculations are simply an extension of the calculations 
of Thurston and Morrison (9) to values of N upto 300. In comparison 
with the earlier calculations of Zimm and Tschoegl, the importance of our 
calculations and those of Thurston and Morrison is twofold: (i) our 
calculations are exact (within the confines of the Zimm treatment of the 
bead/spring model); and (ii) N, the number of springs per polymer molecule, 
appears as an additional parameter whose value is at our disposal when 
we wish to compare theoretical predictions with experimental! data on any 
given solution of linear polymer molecules of equal molecular weight. 

It is possible that (ii) is more important, and more open to question, 
than (i). Thurston and Morrison (9) have already made use of (ii) in offering 
a possible interpretation of the dependence of [n] on M, using the 
assumption that N is proportional to M fora given polymer/solvent system. 
Thurston and Schrag (15) have similarly interpreted their experimental data on 
oscillatory flow birefringence; for values of M down to 10° , they use values 
of N down to I for the polystyrene/Aroclors system. 

We do not propose here to embark on a comprehensive discussion of 


the issues raised by the use of N as a parameter, but we do recognize that 
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this use is open to the criticism that N_ is too closely connected with the 
main artificial feature of the bead/spring model (namely, the arbitrary 
concentration of hydrodynamic interaction at discrete points on the Gaussian 
spring system) to provide a meaningful parameter having significance for 

any actual polymer solution. Awareness of this criticism no doubt influenced 
the early applications of the bead/spring theory (1, 2, 4, 5) in which the 
equations were cast in a form in which N did not appear: this was achieved 
partly taking N =© in the upper limits of sums in equations such as { 3.2], 

[ 3.4], and { 3.6] and partly by the use of the integro-differential eigenvalue 
equation for the case of very large N. On the other hand, the use of very 
large values of N is itself open to criticism: for a given polymer molecule, 

a certain minimum number (perhaps not well-defined, but nonetheless real) of 
backbone bonds are required in order that the corresponding part of the molecule 
may be approximately treated as a Gaussian spring; furthermore, the fact that 
N may not appear in an equation does not, by itself, eliminate the dependence 
of the equation on the artificial feature of the model. 

We are perhaps on safer ground if we confine applications of the bead/ 
spring theory to properties for which the results predicted by the theory are 
insensitive to the value of N; an example is furnished by the low-frequency 
behavior of the dynamic functions shown in Figure 7, On the other hand, it 
is just conceivable that a more fundamental treatment of polymer/solvent 
dynamics than that represented by the bead/spring idealization could furnish 


criteria for the validity of the bead/spring idealization (with finite N) 
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as an approximation. We hope that the calculation of the present report, 


when compared critically with suitable experimental data, might lead to a 
better understanding of the justification for and limitations of the bead/ 


spring model. 
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Table 3: Exact characteristic values and steady-and oscillatory- 
shear-flow functions 


Key to entries 


ETAIR ETA2R GIR 


[4.3] [4.4] uw 
Exponential notation: .769-03 means 0.769 x 10° 


Parameter values 


0.05 0.075 0.1 0.13 0.15 0.2 
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NX 
2 +2: 
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l op 
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